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We propose another integrate-and-fire model as a single neuron model. We 
study a globally coupled noisy integrate-and-fire model with inhibitory interac- 
tion using the Fokker-Planck equation and the Langevin equation, and find 
a reentrant transition of oscillatory states. Intermittent time evolutions of 
I neuron firing are found in strongly inhibited systems. We propose another 

' integrate-and-fire-or-burst model including the dynamics of the low-threshold 

I Ca^+ current based on the new integrate-and-fire model. We study a globally 

' coupled noisy integrate-and-fire-or-burst model with inhibitory interaction us- 

I ing the Fokker-Planck equation, and find bistability of the tonic mode and burst 

1^ ' mode. Doubly periodic oscillation appears in a coupled system of two neuron 

' assemblies, which is similar to the spindle oscillation in thalamic cells. 

o ■ 

^ ! 1 Introduction 

cr. 

^ I The limit cycle oscillation and the synchronization are found in various natural 

• , phenomema such as the firing of fireflies and the stick-slip motion of a fault in the 

' earthquake. Coherent oscillation appears as the global synchronization of the 

^ . coupled oscillators. The synchronization of oscillators is also considered to play 

important roles in neural information processing 0. Winfree and Kuramoto 
studied the global synchronization in general coupled oscillators including the 
phase oscillators The leaky- integrate- and-flre model (IF model) is one 

of the simplest models for a single neuron and often used to study dynamical 
behaviors of neural networks. Each neuron receives an input via synaptic con- 
nections from other neurons. The neuron fires, when the membrane potential 
goes over a threshold, and sends out an impulse as an output to other neurons. 
MiroUo and Strogatz studied a globally coupled system of integrate-and-fire neu- 
rons, and showed that perfect synchronization occurs in a finite time jll. Mutual 
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synchronization is often observed among a large number of neurons with exci- 
tatory interaction. Inhibitory neurons can exhibit synchronization under some 
appropriate conditions However, the synchronization among inhibitory 

neurons is not so robust against some noises or heterogeneities, as was studied 
in several models in Refs. [8] and [9]. On the other hand, the synchronization 
among thalamic neurons is known to be rather robust experimentally, although 
the synaptic coupling is inhibitory. The synchronous firing among the thalamic 
neurons generates characteristic brain waves called spindle waves in an early 
stage of sleep 10 . The reason of the robust synchronization among inhibitory 
neurons is not well understood. In this paper, we will study the synchronization 
among noisy inhibitory neurons, using a new type of IF model. 

A single neuron may be affected by environmental noises, then, the model 
equation is generalized to a Langevin type equation such as a noisy integrate- 
and-fire model. In the noisy system, the coherent oscillation appears as an ana- 
logue of the phase transition in the statistical mechanics. Nykamp-Tranchina 
and Brunei-Hakim showed that a population density approach using the Fokker- 
Planck equation is a useful method to study a large neural network of noisy 
integrate-and-fire neurons ^2 ^1 ^] ■ We studied a nonlocally coupled noisy 
integrate-and-fire model and found a traveling pulse state using the Fokker- 
Planck equation corresponding to the Langevin equation |14| . The phase tran- 
sition in a large population of noisy elements can be studied with the Fokker- 
Planck equation in a clear manner, since the Fokkcr-Planck equation is a deter- 
ministic equation. In this paper, we apply the method to a large population of 
noisy integrate-and-fire neurons to clarify several phase transitions in some mod- 
els of neural networks, although related subjects have been intensively studied 
with different methods and models [T^ITHlirT] . 

The sleep spindle oscillations are a kind of brain waves, which appear in an 
early stage of sleep. The oscillations have a characteristic waveform of waxing- 
and-waning. It is known that these oscillations are generated by the interaction 
between thalamic reticular (RE) and thalamocortical (TC) cells ^1 118L I19| . 
The interactions from RE cells to TC cells are inhibitory and the interactions 
from TC cells to RE cells are excitatory, the mutual interaction among RE 
cells is inhibitory, and there is no mutual interaction among TC cells. The low- 
threshold Ca^"*" current plays an important role for the excitation of the thalamic 
cells. The integrate-and-fire-or-burst model (IFB model) is a generalized model 
of the IF model, which takes the dynamics of the low-threshold Ca^+ current 
into consideration jl8j . The dynamics of the IFB model under external stimuli 
was studied in Ref. [19]. We will propose a new type of IFB model including 
the dynamics for the low-threshold Ca^+ current based on the new IF model. 
We will study a globally coupled noisy IFB model with inhibitory interaction 
as a neural network model of the RE cells, and finally study a coupled systems 
of two types of IFB neuron assemblies as a coupled model of the RE cells and 
TC cells. 

We will investigate the coupled noisy integrate-and-fire models stepwise. In 
§2, we propose a new type of integrate-and-fire (IF) model with inhibitory cou- 
pling and show the synchronization between two neurons. In §3, we investigate 
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a globally coupled noisy integrate-and-fire model using the Fokker-Planck equa- 
tion. In §4, we propose a new type of integrate-and-fire-or-burst (IFB) model 
by combining the new type of integrate-and-fire model with the dynamics of the 
low-threshold Ca^"'' current, and investigate the global oscillation in the globally 
coupled noisy IFB model. In §5, we investigate a coupled system of two neuron 
assemblies and show doubly periodic oscillation. 

2 Synchronization of two integrate-and-fire neu- 
rons with inhibitory coupUng 

In the usual integrate-and-fire model, the membrane potential v is reset to a cer- 
tain value Vr ~ —50 mV, just after the membrane potential reaches a threshold 
value, which is assumed to be about —35 ^ —45 mV. In more realistic models, 
the neuron is excited after the membrane potential goes over the threshold, the 
membrane potential takes positive values for a while, and then the membrane 
potential returns to a negative value around the resting potential. Taking the 
dynamics in the excited state into consideration, we propose another type of 
simple integrate-and-fire model. Our new IF model keeps the most favorable 
point of the IF model that the firing of a neuron can be described using only one 
variable, in contrast to the Hodgkin-Huxley model including four variables |2()j 
or the FitzHugh-Nagumo model including two variables |21|. 

We study the synchronization of two neurons with inhibitory coupling, using 
the new IF model. The model equation is written as 

a{vi - Vo) + Io + Is, for < Vr, 

- Vq) + Is, for Vi>VT, (1) 

where Vi denotes the membrane potential for the ith cell, Vb = — 65 denotes 
the resting potential, Vq is another parameter which controls the dynamics in 
the excited state Vi > Vt and is fixed to be Vq = 35 in this paper, Iq is a 
constant external stimulus current. Is is a synaptic current, and C = 2 and 
a = 0.035 are constants. The units of v and t are mV and ms. When the 
membrane potential exceeds the threshold Vr — —35, the cell is excited. In the 
usual integrate-and-fire model, Vi is reset to a potential Vr — —50 just after 
the firing. However, in our generalized model, the membrane potential jumps 
to a potential Vi > (Vi = 60 is assumed in this paper.) and then it tends to 
decay toward Vq if Is = 0, obeying the second equation of Eq. (1). When the 
membrane potential reaches a second threshold V2, Vi is reset to Vr. The two 
kinds of time evolutions correspond to the slow dynamics on the two branches 
of nuUclines in the excitable equations such as FitzHugh-Nagumo model [?T] 
and the McKean model |22| . 

There are several model equations which describe the dynamics of the synap- 
tic current Is{t) |^3,. In this paper, we use a simple model including a response 
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Figure 1: (a) Two time evolutions by Eq. (1) starting from vi = —50 (solid 
line), V2 = 50 (dashed line) and h = -0.9973 for Iq = 1.5, K = -0.5, V2 = 40 
and T = 10. (b) Stability exponent r as a function of t ioi K = 0.1, —0.1 and 
—0.5 at Iq = 1.5 and r — 10. (c) Two time evolutions by Eq. (1) starting 
from vi — —48 (solid line) and V2 — 59.9 (dashed line) and Is — —0.9973 for 
Iq = 1.5, K ~ —0.5, V2 — 58 and r — 10. (d) Stability exponent r as a function 
of V2 for K = 0.1, -0.1 and -0.5 at /q = 1.5 and r = 10. 

time. The synaptic current for a coupled system of two neurons is assumed to 
obey 

T'^^-{Is-K/2j2e{v,)). (2) 

where r is a time constant, K < Q in an inhibitory coupling constant, and 0{x) 
is the Heaviside step function. If the cell is excited {vj > 0), the inhibitory 
synaptic current Is < is induced. This type of simple synaptic coupling is 
often used in some theoretical models j24) . 

Figure 1 (a) displays the time evolution of vi and V2 starting from vi (0) = 
-50, V2{0) = 50 and 7,(0) ^ -0.9973 at Iq = 1.5, K = -0.5, V2 = 40 and 
T = 10. The two oscillatory time evolutions approach each other, which means 
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that the mutual synchronization is attained. The solution of the synchronized 
oscillation can be explicitly solved, because the model equation is a piecewise 
linear equation. The stability of the synchronized motion can be studied by the 
linearized equation of Eqs. (1) and (2). The deviation Sv = fi — V2 obeys a 
simple equation: Cd5v/dt = —a5v. The deviation 5v changes discontinuously 
at the moment t\ of the first threshold satisfying v\{t{) = Vt, as dv{ti+) = 
{dv / dt)^-^^Vi / {dv / dt)^-^^VT^^{^i-) 1 ^-iid at the moment t2 of the second thresh- 
old satisfying V\{t2) = V2, as Sv{t2+) = {dv/dt)y^=VR/{dv/dt)y^=v2^v{t2-)- 
Therefore, the deviation Sv{t) grows Sv{t + T) = r6v{t) after a period T of 
the synchronized motion, where the stability exponent r is expressed as 

{dv/dt)y,=v, idv/dt)y,^v^ ^_^^c!-T 

{dv/dt)y^=VT {dv/dt)y^=v2 

Figure 1(b) displays the stability exponent r as a function of r for Iq = 1.5, K = 
—0.5 and r = 10. If the stability exponent r is smaller than 1, the synchronized 
motion is stable. As K <0 and \K\ is larger, the synchronized motion becomes 
more stable. There is an optimum response time t where r takes a minimum 
value, although the synchronized motion is stable even for very small r for 
K < 0. In the case of excitatory coupling ii' > 0, r is larger than 1 and the 
synchronized motion is unstable. That is, the mutual synchronization is attained 
in the case of inhibitory coupling and the synchronized motion becomes more 
stable, as the inhibitory coupling is stronger. 

If the second threshold V2 approaches Vi, the duration time of the excita- 
tion becomes smaller, and in the limit of V2 = Vi, the duration time of the 
excitation becomes and our model approaches the original integrate-and-fire 
model. Figure 1(c) displays the mutual synchronization of two neurons start- 
ing from vi = —48 (solid line) and V2 = 59.9 (dashed line) and = —0.9973 
at Iq = 1-5, K = —0.5, V2 = 58 and r = 10. The duration time of the ex- 
cited state is much shorter than the case of Fig. 1(a) owing to the smallness 
of |Vi — V2I, however, the mutual synchronization is still attained. Figure 1(d) 
displays the stability exponent r of the synchronized state as a function of V2 
for K = —0.5; —0.1 and 0.1 at Iq = 1.5 and r = 10. As V2 approaches Vi = 60, 
the stability exponent approaches 1, because the duration time of the excited 
state becomes shorter, and then the synaptic interaction becomes weaker ef- 
fectively. In the following numerical simulations, we use the parameter value 
of V2 = 40 (except in Fig. 4(b)), because mutual synchronization occurs fairly 
strongly, although the duration time of the excitation might be too long as a 
typical neuron. 

3 Global synchronization in a large population 
of noisy integrate-and-fire neurons 

Next, we consider a large population of integrate-and-fire neurons with in- 
hibitory coupling. The total number of the neurons is assumed to be A^. The 
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elemental dynamics for each neuron obeys Eq. (1). The synaptic current is 

expressed as 

rIT ^ 
i=i 

A completely synchronized state where Vi{t) = v{t) for any i is a special solution 
to the globally coupled system, and the stability exponent r of the completely 
synchronized state is calculated from the linear growth rate of 6vi (t) = Vi {t) — 
v{t). The stability exponent r takes the same value as that in a coupled system 
of two neurons, because K/NJ2j'=i (^{'"j) and Is(t) take the same value as the 
case of two neurons. That is, the completely synchronized state is stable for the 
globally coupled system with inhibitory interaction without noise terms. 

If a noise term is added to Eq. (1), the model equation is expressed as a 
Langevin equation: 

(11) ' 

= -a{vi - Vo) + C^t), for Vi < Vt, 

dv ' 

C-^ = -a{v, - V^) +h + C^i{t), for Vi > Vt, 

where^i(t) denotes Gaussian white noise which satisfies {Ci(t)£,j(f/)) = 2DSij5{t— 
t'). The corresponding Fokker-Planck equation for ^ cx) is written as 

^ = -^{f{xJoJs)P{x)] + D^+5{x-Vr)Ji{t) + 6{x-VR)J2{t), 
= -{Is-K P{x)dx), (5) 

where P{x, t) denotes the probability density that the membrane potential v 
takes a value a;, Ji{t) = —D{dP/dx)x=VT and J2{t) — —D{dP/dx)x=v2- The 
function f[x,Io,Is) denotes a simplified form of the deterministic part of the 
elemental dynamics for each neuron, that is, the first two equations in Eq. (4) 
are rewritten as dvi/dt = f{vi, /q, Is)+S,i{f^ for the sake of simplicity. The terms 
including the (5-function represent the two jmnp processes from w = Vr to Vi 
and from V2 to Vr. The probability density P{x) is at a; = Vt and V2. The 
integral /g P(x)dx denotes the fraction of the excited neiirons, which is equal 

to Yl^j^i^i'^j)/^ ■ We have performed direct numerical simulation of Eq. (5) 
with the Euler method with Aa; = 0.1 and At = 0.005. Similar direct numerical 
simulation was performed for the usual integrate-and-fire neurons with only one 
threshold in Ref. [14]. 

We have calculated the average value of the membrane potential E = P{x, t)xdx. 
Figure 2(a) displays time evolutions of E{t) at K = —2 (solid line) and —10 
(dashed line) for /q = 2.5, D = 0.2 and t = 20. Clear limit cycle oscillations of 
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Figure 2: (a) Time evolution of the average membrane potential E{t) by Eq. (5) 
at if = -2 (solid line) and K = -10 (dashed line) for Iq = 2.5, D = 0.2, V2 = 40 
and T ~ 20. (b) Time evolution of the probability density P{x, t) a,t K — —2. 
Since P{x, t) is always in —30 < x < 30, the region is not plotted, (c) Time 
evolution of the probability density P{x,t) aX K = —10. 



the average membrane potential are observed, which implies the synchronization 
of firings in noisy environment. At K = —2, the average membrane potential 
goes up to i? ~ 35, but at if = —10, the average membrane potential is always 
below 0. Figures 2(b) and (c) display time evolutions of the probability density 
P{x,t) at (b) K = -2 and (c) K = -10 for Iq = 2.5, D = 0.2 and r = 20. 
The probability density is oscillating across the threshold a; = Vr as a whole 
at K = —2, but only a fraction of probability density goes up to the positive 
region x > at K = —10, remaining most part of probability density below the 
threshold. That is, almost all neurons fire synchronously at if = —2, but only a 
fraction of neurons fire synchronously at if = —10, which will be clearly shown 
by the corresponding Langevin simulation in Fig. 4. 

Figure 3(a) displays the peak-peak amplitude AE — Emax ~ Emm of the 
oscillation of the average membrane potential E{t) for Iq = 2.5, t = 20 and 
D = 0.2 as a function of the coupling constant |if | = —if. As is seen in 
Fig. 3(a), the global synchronization occurs for 1.1 < |if | < 3.4 and |if | > 8.3. 
The transitions are supercritical Hopf bifurcations. The global synchronization 
reappears in a parameter region of smaller |if |. That is, a reentrant transition 
of the global synchronization occurs in this noisy system. Figure 3(b) shows a 
parameter region where the global synchronization is observed in a parameter 
space of r and |if | for Iq = 2.5 and D = 0.2. There are two parameter regions 
where the global oscillation occurs. In the parameter region of smaller |if|, 
almost all neurons fire synchronously. The parameter region of the smaller |if | 
is bounded between 6 < t < 40. This might be related to a fact that the stability 
exponent r is the smallest near r = 15 for Iq — 1.5 as shown in Fig. 1(b), that is, 
the finite response time facilitates the global synchronization. In the parameter 
region of larger |if |, a fraction of neurons fire synchronously, but the average 
membrane potential is lower than 0. The strong inhibition with large if inhibits 
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Figure 3: (a) Oscillation amplitude as a function of \K\ for /q = 2.5, D = 0.2 
and r = 20. (b) Parameter region where the global synchronization appears in 
the parameter space of (r, \K\) for Iq = 2.5 and D = 0.2. The coupling constant 
K is negative, and the same parameter values Iq and D are used for (a) and (b) 



the firing of all neurons. Neurons which fire earlier, inhibit the firing of the other 
neurons. Which neuron fires is not determined but change randomly in our noisy 
system, in contrast to a clustered state found in some deterministic systems |25) . 
That is, the global oscillation appears, but the globally synchronous firing does 
not occur in the region of large K. 

We have shown numerical results for the Fokker-Planck equation. The 
Fokker-Planck equation is favorable to investigate the global oscillation clearly, 
since the equation is deterministic. However, the time evolution of individual 
neuron cannot be seen. The direct Langevin simulation can show the time evo- 
lutions of each neuron. The difference of the two types of global oscillation can 
be clearly seen in the Langevin simulation of Eq. (4). We have performed a 
direct Langevin simulation for N — 2000 neurons, and got time evolutions of 
membrane potentials w^'s for all neurons and the average value of the membrane 
potentials. We have performed some numerical simulations by changing parame- 
ter values of D, V2 and K. Figures 4(a), (b) and (c) display the time evolutions 
of Vi^i{t) of the first neuron (solid line) and the average {v) = "^jLiVj /N 
(dashed fine) for Iq = 2.5, r 20 and N — 2000. Figure 4(a) displays the time 
evolutions at K — —2, D = 0.2 and V2 — 40, which corresponds to the solid line 
in Fig. 2(a). The firing of each neuron is almost synchronized to the average 
activity (v) . The timing of the firing is slightly different for each neuron owing 
to the external noise, that is, the phase of the oscillator is randomized by the 
external noise. The average value (v) exhibits limit cycle oscillation, which is al- 
most equal to E{t) in Fig. 2(a). The global synchronization may be interpreted 
by the concept of the phase synchronization as studied in Ref. [3]. Figure 4(b) 
displays the time evolution at K = —2, D — 0.02 and V2 = 58. At V2 = 58, the 
duration time of the excited state is much shorter, and the effective interaction 
becomes weak as shown in Fig. 1(c) and (d), but almost synchronized firing is 
seen for the smaller noise strength D = 0.02. For the same parameter values 
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Figure 4: Time evolutions of Vi (solid line) and the average (v) (dashed line) 
in the Langevin simulation of Eq. (4) for a globally coupled noisy IF model 
at (a) K ^ -2,D ^ 0.2, V2 = 40, (b) K = -2, D = 0.02, V2 = 58 and (c) 
K = -10, £1 = 0.2,1/2 ^ 40 for Iq ^ 2.5, t ^ 20 and N = 2000. (d) Raster 
plot at lo = 2.5, K = -10, D = 0.2, V2 40,r = 20 and N = 2000 which 
corresponds to (c) , in which a dot is plotted at time t = rt (rt is an integer) and 
the position i satisfying Vi{t = n) > for 1 < i < 100. 



K — —2 and D = 0.2 as Fig. 4(a), the global synchronization was not observed 
at V2 = 58, because the timing of the firing is randomly distributed. Inhibitory 
neurons can exhibit the global synchronization even if the duration time of the 
excited state is short, but the stability becomes weaker. Figure 4(c) displays 
the time evolutions at K = —10,D — 0.2 and V2 — 40, which corresponds to the 
dashed line in Fig. 2(a). The average value (v) exhibits limit cycle oscillation, 
which is close to E{t) in Fig. 2(a). The activity Vi{t) of each neuron shows that 
the firing occurs intermittently at some peaks of {v). A fraction of neurons fire 
near the peak of (v), and the other neurons are inhibited to fire by the strong 
global inhibition by the fired neurons. Which neurons are selected for firing is 
not determined but changes randomly. As a result, the intermittent time evo- 
lution appears as shown in Fig. 4(c). The random but somewhat synchronized 
firing is seen in the raster plot of Fig. 4(d) at the same parameter values as in 
Fig. 4(c), where firing sequences of neurons for 1 < i < 100 are denoted by dot 
patterns. Intermittent firing was found also in a numerical simulation by Wang 
et al.|16|. but their model is deterministic and the coupling strengths among 
neurons are randomly distributed. Similar intermittent behaviors are observed 
also experimentally |2(>j . 
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4 Bistability in a large population of the integrate- 
and-fire-or-burst neurons 



The integrate-and-fire model is a simpler model than realistic model equations 
such as the Hodgkin-Huxley equation 20 , in that the detailed dynamics of ion 
channels are neglected. Specific ion channels play an important role in the fir- 
ing process of specific neurons. The low-threshold Ga?^ current is important 
for thalamic neurons. The integrate-and-fire-or-burst model was proposed as 
a simple model to describe the dynamics of the membrane potential of thala- 
mic cells 18 . It is a coupled equation of the integrate-and-fire model and the 
dynamics of the low-threshold Ca^+ current Ixit)- The current Irit) is given 
by Irit) — gTh{t){VH — '^{t))S{v{t) — Vh), where the parameters are given as 
qt — 0.07, Vh — 120 and Vh = —60, and h{t) is a slow variable determined by 
the membrane potential v. In the integrate-and-fire-or-burst model, the slow 
variable h{t) (0 < ft. < 1) is assumed to obey 

Mv)^ = -h + eiVh-v{t)), 

where Th = 20 for v > Vh and Th = 100 for u < 14 [1111111. The low-threshold 
current It is zero, if the membrane potential is constant. The current It fiows 
in the neuron for a short time Tfi, after the membrane potential v goes over 
the threshold Vh from below. The integrate-and-fire-or-burst model is a coupled 
equation of the standard integrate-and-fire model and the equation of h(t). We 
study a new integrate-and-fire-or-burst (IFB) model, which is a coupled equation 
of the new IF model Eq. (1) and the equation of h{t). We investigate a globally 
coupled noisy IFB model, which is expressed as 

C-^ = -a{v,-VQ) + I^,+Is + IT^ + CUt), ior v,<Vt, 
du ■ 

= -a{v, - V^) + /, + IT^ + C^,{t), for > Vt, 

Th^ = -K + 9iVh-vS)), lT^{t)^gM){VH~vS))e{v,{t)-Vh), 
dt 

The slow variable hi{t) is another stochastic variable, because Vi{t) is a 
stochastic variable. Therefore, the corresponding Fokker-Planck equation takes 
a form of a two-dimensional partial differential equation: 

^ ^ -^if('^'^o,Is+lT)}P{x)-^[{{-y + 9{Vh-x)/rh{x)}P] 
+ + '^(^ - Vi)Mt, y) + Six - VR)J2{t, y), 

- -i^s-Kj^ j^P{x)dxdy), (7) 
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where P(x,y,t) is the probability density that the membrane potential v takes 
a value x and the slow variable h takes a value y, f{x,Io,Is + It) denotes 
the deterministic part of the elemental dynamics of the membrane potential, 
lT{x,y) = gTyiVH-x)0{x-Vh), and Ji{t,y) = -D{dP/dx)x=VT and J2 (*,?;) = 
—D{dP/dx)x=V2- Since h{t) does not receive a random force, there is no diffu- 
sion term with respect to y in Eq. (7). The deterministic forces to the membrane 
potential v{t) and the variable h(t) appear in the first two drift terms on the 
right-hand side of Eq. (7). We have performed numerical simulations of the 
two-dimensional partial differential equation (7) using the Euler method with 
Ax = 0.1, At = 0.005 and Ay = 0.05. 




Figure 5: (a) Two time evolutions of the average membrane potential E{t) by 
Eq. (7) for Iq = 1.6,K = -40, D = 0.2 and r = 20. The initial conditions for 
the two time evolutions are different, (b) The range of the temporal variation 
of the average value: (/i) as a function of Iq for K = —40, D = 0.2 and r = 20. 
The solid line denotes the burst mode and the dashed line denotes the tonic 
mode, (c) Oscillation amplitude of E{t) as a function of Iq for the burst mode 
at K = —40 and D = 0.2. (d) Oscillation amplitude of E{t) as a function of Iq 
for the tonic mode at K = —40 and D = 0.2. 

Figure 5(a) displays two time evolutions of E(t) = P{x, y)xdxdy for 

Iq = 1.6, K = —40,1? = 0.2 and r = 20. The initial conditions are different 
for the two time evolutions. That is, the two kinds of stable oscillatory states 
have appeared for the same parameters. It implies the bistability. The average 
membrane potential is always larger than Vh in the dashed line, so h{t) is almost 
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0. The low-threshold Ca^+ current docs not play a role, and the frequency of the 
oscillation is relatively high. This is called the tonic mode. On the other hand, 
the average membrane potential E{t) oscillates around Vh = —60 in the solid 
line. The low-threshold Ca^+ current plays an important role in this mode. The 
low-threshold current It depolarizes the IFB neuron and the firing is facilitated. 
If many neurons are fired, the strong inhibitory coupling hyperpolarizes each 
nciiron and the membrane potential goes down below Vh, and then is rebound to 
go over the threshold again. The low-threshold Ca^+ current facilitates the 
oscillation of the firing, and therefore the amplitude of oscillation is relatively 
large and the frequency of oscillation is relatively low in this mode. This mode is 
called the burst mode or the rebound mode. The tonic mode and the burst mode 
are bistable at the parameters in Fig. 5(a) in our model. Figure 5(b) displays the 
range of the temporal oscillation of the average {h{t)) ~ Jq Pi^, y, t)ydxdy 
as a function of Jq for K = —40, D = 0.2 and t = 20. The finite width of the 
temporal variation of {h{t)) at a certain Iq implies the limit cycle oscillation 
of {h{t)). The upper branch represents the burst mode and the lower branch 
{{h) = 0) represents the tonic mode. (We have distinguished the two modes 
in this paper from a viewpoint that the low-threshold Ca^+ current plays an 
important role or not, that is, the average value of h{t) is finite or almost zero.) 
The bistability occurs for 0.55 < /q < 2.3. Figure 5(c) displays the peak-peak 
amplitude /S.E of the burst oscillation as a function of 7o- The supercritical 
Hopf bifurcation is seen at /q ^ —0.15. Near the threshold, the oscillation 
is sinusoidal, but a period doubling bifurcation occurs for the burst mode for 
Iq > 0.7. The burst-mode oscillation shown by the solid line in Fig. 5(a) displays 
the time evolution of E{t) after the period doubling bifurcation. In the tonic 
mode, the Hopf bifurcation occurs near Jq = 0.95, as shown in Fig. 5(d). Figure 
5(c) and (d) show that the oscillation amplitude in the burst mode is relatively 
large compared to the tonic mode. 

5 Spindle oscillation in a coupled system of two 
neuron assemblies 

We investigate a coupled system of two neuron assemblies in this section, as a 
model of the thalamic neural network. There are two kinds of neurons: thalamic 
reticular (RE) cells and thalamocortical (TC) cells in the thalamus. The low 
threshold Ca^+ current plays an important role in both types of thalamic cells. 
The interaction among RE cells is inhibitory. In the previous section, we have 
investigated the dynamical behaviors of a globally coupled system of one kind of 
inhibitory IFB neurons, which corresponds to an assembly of RE cells alone. It is 
known that the interaction between the RE cells and the TC cells is important 
for generation of spindle waves, which are characteristic brain waves in the 
second stage of sleep. We will study a coupled system of two neuron assemblies 
corresponding to RE cells and TC cells. It is known that the interaction from 
the RE cells to the TC cells is inhibitory and the interaction from the TC cells 
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to the RE cells is excitatory, and there is no mutual interaction among the TC 
cells. We investigate a Fokker-Planck equation for the coupled system of two 
neuron assemblies composed of noisy IFB neurons. The Fokker-Planck equation 
for the coupled neuron assemblies is expressed as 

BP 3 3 

^ -Q^{fi^^l0lJsl+Is2+lT)Pl}-Q^[{{-y + e{Vhl-x))/Th}Pl] 

+ + 5{x - Fi)Jii(<, y) + 5{x - Vn)Ji2{t,v), 

3P 3 3 

-gf = -g^{fix,l02,Is3+lT)P2}-g^[{i-y + e{V,2~x))/r,}P,}] 
.d^P2 



disi 



n 


dt 




dls2 


T2 


dt 




dls3 


T"3 


dt 



+ '^(^ - yi)J2iit^ y) + - yR)J22{t, y), 

-{Isi-Ki / / Pi{x)dxdy), 

Jo Jq 

-{Is2-K2 / / P2{x)dxdy), 

Jo Jo 

pOO pi 

~{Is3~K3 / / Pi{x)dxdy), (8) 



Jo Jo 

where Pi (x, y, t) and P2 {x, y, t) are the probability densities for the two neuron 
assemblies, Jii{t,y) = -D{3Pi/dx)^^VT ,Ji2{t,y) = -P'{dPi/dx)^^v2,J2i{t,y) = 
—D{dP2/dx)x=VT and J22(i, y) — —D{dP2/ dx)x=v2- Our model is a mathemat- 
ical model and may not describe the dynamics of realistic thalamic cells well, 
but we call the first neuron assembly the RE neurons and the second neuron 
assembly the TC neurons, for the sake of simplicity. We assume the sign of the 
three coupling constants as Ki < 0, K2 > and K3 < 0, taking the thalamic 
neural network into consideration. The synaptic current Isi < represents the 
inhibitory coupling among the RE neurons, Is2 > represents the excitatory 
coupling from the TC neurons to the RE neurons, and Ig^ < represents the 
inhibitory coupling from the RE neurons to the TC neurons. The external in- 
puts to the RE and TC neurons are respectively denoted by Iqi and Iq2. The 
response time Ti for the RE neurons is assumed to be smaller than the response 
times T2 and for the interaction between the RE and TC neurons, since the 
distance between the RE neurons and TC neurons is larger than the distance 
among the RE neurons. The threshold values of the low-threshold Ca^+ cur- 
rent for the RE neuron and the TC neuron are respectively assumed to be 
Vhi = —60 and Vh2 = —70 ^Hl- This model equation is rather complicated 
and includes many parameters. We have performed numerical simulations of 
Eq. (8) for n = 5, T2 = = 100, Ki ^ K3 ^ -30, K2 = 30, £> = 0.2 as an 
example. Figure 6 displays the time evolutions of the average membrane poten- 
tial Ei{t) = Pi{x, y, t), xdxdy of the RE neurons for (a) Iqi = I02 = 2.2 
and (b) Iqi = —0.25, /02 = 1, and the average membrane potential E2{t) of the 
TC neurons for for (c) Iqi = I02 = 2.2 and (d) Iqi = -0.25, /02 = 1- Doubly- 
periodic oscillations are seen in the time evolutions of Ei{t). The long-term 
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Figure 6: Time evolutions of the average membrane potential Ei{t) of the RE 
neurons in Eq. (8) at (a) Iqi = /02 = 2.2 and (b)/oi = —0.25, /02 = 1, and the 
corresponding average membrane potential E2(t) of the TC neurons at (c)/oi = 
I02 = 2.2 and (d)/oi = -0.25, I02 = l- 



oscillatory motion is caused by the mutual interaction between the RE and TC 
neurons. The fast oscillation appears owing to the synchronization among the 
RE neurons, which was shown in the previous section. The two kinds of os- 
cillations are overlapped and the doubly-periodic motion appears in the time 
evolution of Ei{t). The mechanism of the doubly periodic motion is as follows. 
The fast oscillation of the RE neurons repeats for a while. The TC neurons are 
inhibited strongly below the Vh2 by the repetition of the fast oscillations of the 
RE neurons. The firing of the RE neurons becomes weak, as the input from 
the TC neurons becomes small. The TC neurons are rebound after the strong 
inhibition, and the membrane potential of the TC neurons goes over Vh2- The 
low threshold Ca^+ current flows into the TC neurons and the TC neurons fire. 
Then, the RE neurons are excited by the firing of the TC neurons. The long 
period is determined by the slow rebound time. The waxing-and- waning motion 
of the fast oscillation is seen clearly in Fig. 6(b). The fast oscillations repeat 
several times and then disappear. In the repetition of the fast oscillations, the 
period increases gradually, the amplitude of the fast oscillation is decreasing. 
This type of motion is characteristic of the spindle oscillation in the thalamic 
networks j7. . 

In our model, the doubly periodic oscillation is caused by the interaction 
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between the RE neurons and TC neurons, but there are other neuron models 
where the slow oscillation is caused by slow dynamics of some ion channels [211 
■ In their models, the doubly periodic behavior can be generated even in one 
neuron. 




Figure 7: Time evolutions of Vi (solid line) and the average (v) (dashed line) 
in the Langevin simulation corresponding to Eq. (8) with ti = 5, T2 ^ = 
100, Ki=K3^ -30, K2 = 30, £» = 0.2, /qi = -0.25, /02 = 1, and TV = 1000. 



The spindle oscillation is also seen in the corresponding Langevin simulation. 
Figure 7 displays the time evolutions of the membrane potential vi of the first 
RE neuron and the average value (v) for the RE neurons at ti — 5, T2 = = 
100, Ki= K3 = -30, K2 = 30, £» = 0.2, /qi = -0.25, /02 = 1 and TV = 1000, 
which corresponds to Fig. 6(b). The time evolution of the average value (v) 
exhibits the doubly periodic oscillation similar to Fig. 6(b). The firing of each 
neuron is well synchronized to the averaged membrane potential, but the firing 
is sometimes skipped and the intermittent firing is seen. The intermittent firing 
might be characteristic of strongly inhibited systems. Similar intermittent firing 
of the RE neurons is observed experimentally ■ 

The bistability is also observed in this coupled system. Figure 8(a) displays 
the maximum value in the time evolution of the average (hi {t)) — Jq Pi {x, y, t)ydxdy 
as a function of /q = /oi = Io2- The upper branch represents the burst mode 
and the lower branch represents the tonic mode for the RE neurons. The bista- 
bility occurs in 0.3 < Iq < 1.4. Figure 8(b) displays the maximum value in the 
time evolution of the average (/ii(t)) — Jg Pi{x^ y, t)ydxdy as a function of 
/oi for a fixed value of /02 = 1- The bistability is observed in 0.5 < /qi < 1.8, 
although the classification to the two burst and tonic modes seems to be difficult 
in this case. The thalamus is considered to be a relay point between the sensory 
organs and the cerebrum. The bistability of the activity of the coupled system 
of the RE cells and TC cells might be related to some switching function of the 
information transmission at the relay point. 
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Figure 8: Maximum value of the average {hi{t)) = Pi{x,y,t)ydxdy as 

a function of (a) Iq = Iqi = Ia2, and (b) as a function of /qi at I02 = 1- Both 
figures display the existence of the bistability and the hysteresis. 



6 Summary 

We have proposed another IF model including the dynamics in the excited state. 
We have studied the stability of the synchronized motion of a coupled system of 
two neurons and N neurons. In our model, the synchronization is rather robust 
even in the case of inhibitory coupling. 

We have investigated global synchronization of noisy IF neurons with in- 
hibitory coupling using the Fokker-Planck equation, and found the reentrant 
Hopf bifurcation. In the parameter region of weak coupling, almost all neurons 
fire synchronously, but in the parameter region of strong coupling, only a frac- 
tion of neurons fires synchronously owing to the strong inhibition. Intermittent 
firing is seen in each neuron as is clearly shown by the Langevin simulation. We 
have proposed another IFB model including the dynamics of the low-threshold 
Ca^+ current based on the new IF model. We have performed direct numer- 
ical simulations of the Fokker-Planck equation for the globally coupled IFB 
model, and found the bistability of the tonic model and the burst mode. In the 
burst mode, the low-threshold Ca^+ current facilitates the global oscillation. 
We have further performed numerical simulations of the Fokker-Planck equa- 
tion for the coupled system of the RE neurons and TC neurons, and found a 
doubly-periodic motion similar to the spindle oscillation observed in the second 
stage of the sleep. The doubly-periodic motions are also bistable in a certain 
parameter range of the inputs, which might be related to a switching function 
of information transmission at the relay point. 

We have found various bifurcations including the Hopf bifurcation, period 
doubling bifurcation and jump phenomena between bistable states in globally 
coupled noisy IF and IFB models. 

Our model includes many control parameters. Our choice of the parameters 
may not be physiologically suitable. In future work, we will adjust the parameter 
values and perform numerical simulations more relevant to realistic neurons. 
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